%Using Bootstrap method with 500 samples to test MA rules
%input: daily returns, number of days of short term, number of days of long term, band
%output: the fractions of the 500 simulations that result in larger statistics than those for real data
function [ARparameters ARparaTStat ARfractions GARparameters GARparaTStat GARfractions EGARparameters EGARparaTStat EGARfractions] = testMABootstrap(...
    daily_price, risk_free, cost, s_term, l_term, band, option, day_minmax, test_type)

    %Check if the MAtest is FMA, VMA or TRB
    %If test_type = 'F', the test is FMA test, if type = 'V', the test is VMA test, if type = 'T', the test is TRB test
    if test_type == 'F'
        MAtest = @testFLMA;
    elseif test_type == 'V'
        MAtest = @testVLMA;
    elseif test_type == 'T'
        MAtest = @testTRB;
    end
        
    %Calculate the statistics: mean buy, mean sell, mean buy-sell difference, standard deviations for buy and sell returns
    [r, rbuy, rsell] = MAtest(daily_price, risk_free, cost, s_term, l_term, band, option, day_minmax);
    rbuy = rbuy((isnan(rbuy)==0));
    rsell = rsell((isnan(rsell)==0));
    S0 = daily_price(1);
    buy = mean(rbuy);
    sell = mean(rsell);
    buy_sell = buy - sell;
    b_sig = std(rbuy);
    s_sig = std(rsell);
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     ari = arima(1,0,0);
%     [ari vc1] = estimate(ari,r,'print',false);
%     ARsim = simulate(ari,length(r),'numPaths',500);
%     ARparameters = [ari.Constant ari.AR];
%     ARparaTStat =  [ari.Constant/sqrt(vc1(1,1)) cell2mat(ari.AR)/sqrt(vc1(2,2))];
%     ARprice = zeros(length(r)+ 1, 500);
%     ARprice(1,:) = S0;
%     for i = 2:size(ARprice,1)
%         ARprice(i,:) = ARprice(i-1,:).*exp(ARsim(i-1,:));        
%     end     

    y = r(2:end);
    % y = [y];
    x = r(1:end-1);
    x = [ones(length(x),1) x];

    [b,bint,residual,rint,stats] = regress(y,x);

    for i=1:500
        ARsim(1,i) = r(1);
        ARprice(1,i) = daily_price(1);
        noise = datasample(residual,length(residual),'Replace',true);
        for j=2:length(r)
            ARsim(j,i) = b(1)+ b(2)*ARsim(j-1,i) + noise(j-1);
            ARprice(j,i) = ARprice(j-1,i).*exp(ARsim(j-1,i));
        end
        ARprice(j+1,i) = ARprice(j,i).*exp(ARsim(j,i));
    end
     
    ARparameters = b;
    ARparaTStat =  [b(1)/sqrt(var(residual)) b(2)/sqrt(var(residual))];


    gar = arima('ARLags',1,'Variance',garch(1,1));
    [gar vc2] = estimate(gar,r,'print',false);
    GARsim = simulate(gar,length(r),'numPaths',500);
    GARparameters = [gar.Constant gar.AR gar.Variance.Constant gar.Variance.GARCH gar.Variance.ARCH];
    GARparaTStat =  [gar.Constant/sqrt(vc2(1,1)) cell2mat(gar.AR)/sqrt(vc2(2,2)) gar.Variance.Constant/sqrt(vc2(3,3)) ...
        cell2mat(gar.Variance.GARCH)/sqrt(vc2(4,4)) cell2mat(gar.Variance.ARCH)/sqrt(vc2(5,5))];
    GARprice = zeros(length(r)+ 1, 500);
    GARprice(1,:) = S0;
    for i = 2:size(GARprice,1)
        GARprice(i,:) = GARprice(i-1,:).*exp(GARsim(i-1,:));        
    end
    
    egar = arima('ARLags',1,'Variance',egarch(1,1));
    [egar vc3] = estimate(egar,r,'print',false);
    EGARsim = simulate(egar,length(r),'numPaths',500);
    EGARparameters = [egar.Constant egar.AR egar.Variance.Constant egar.Variance.GARCH egar.Variance.ARCH egar.Variance.Leverage];
    EGARparaTStat =  [egar.Constant/sqrt(vc3(1,1)) cell2mat(egar.AR)/sqrt(vc3(2,2)) egar.Variance.Constant/sqrt(vc3(3,3)) ...
        cell2mat(egar.Variance.GARCH)/sqrt(vc3(4,4)) cell2mat(egar.Variance.ARCH)/sqrt(vc3(5,5)) cell2mat(egar.Variance.Leverage)/sqrt(vc3(6,6))];
    EGARprice = zeros(length(r)+ 1, 500);
    EGARprice(1,:) = S0;
    for i = 2:size(EGARprice,1)
        EGARprice(i,:) = EGARprice(i-1,:).*exp(EGARsim(i-1,:));        
    end
    
    count = zeros(3,5);
    
    %count the number of simulations which have the statistics larger than real data's
    for i = 1:500
       [~, s_rbuy1,s_rsell1] = MAtest(ARprice(:,i), risk_free, cost, s_term, l_term, band, option, day_minmax);
       [~, s_rbuy2,s_rsell2] = MAtest(GARprice(:,i), risk_free, cost, s_term, l_term, band, option, day_minmax);
       [~, s_rbuy3,s_rsell3] = MAtest(EGARprice(:,i), risk_free, cost, s_term, l_term, band, option, day_minmax);
       
       s_rbuy1 = s_rbuy1((isnan(s_rbuy1)==0));
       s_rsell1 = s_rsell1((isnan(s_rsell1)==0));
       s_rbuy2 = s_rbuy2((isnan(s_rbuy2)==0));
       s_rsell2 = s_rsell2((isnan(s_rsell2)==0));
       s_rbuy3 = s_rbuy3((isnan(s_rbuy3)==0));
       s_rsell3 = s_rsell3((isnan(s_rsell3)==0)); 
       
       calculations = [mean(s_rbuy1) mean(s_rsell1) mean(s_rbuy1)-mean(s_rsell1) std(s_rbuy1) std(s_rsell1);...
                       mean(s_rbuy2) mean(s_rsell2) mean(s_rbuy2)-mean(s_rsell2) std(s_rbuy2) std(s_rsell2);...
                       mean(s_rbuy3) mean(s_rsell3) mean(s_rbuy3)-mean(s_rsell3) std(s_rbuy3) std(s_rsell3)];
       
       count(:,1) = count(:,1) + (calculations(:,1) > buy);
       count(:,2) = count(:,2) + (calculations(:,2) > sell);
       count(:,3) = count(:,3) + (calculations(:,3) > buy_sell);
       count(:,4) = count(:,4) + (calculations(:,4) > b_sig);
       count(:,5) = count(:,5) + (calculations(:,5) > s_sig);
               
    end
    
    %calculate the fractions
    fractions = count./500; 
    ARfractions = fractions(1,:);
    GARfractions = fractions(2,:);
    EGARfractions = fractions(3,:);
    
        
    ARparameters = ARparameters;
    GARparameters = cell2mat(GARparameters);
    EGARparameters = cell2mat(EGARparameters);
end